{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Use Gustafsson’s Method to Reduce Edge Artifacts](https://swharden.com/blog/2020-09-23-signal-filtering-in-python/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Low-Pass Filter\n",
    "import numpy as np\n",
    "import scipy.signal\n",
    "import scipy.io.wavfile\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def lowpass(data: np.ndarray, cutoff: float, sample_rate: float, poles: int = 5):\n",
    "    sos = scipy.signal.butter(poles, cutoff, 'lowpass', fs=sample_rate, output='sos')\n",
    "    filtered_data = scipy.signal.sosfiltfilt(sos, data)\n",
    "    return filtered_data\n",
    "\n",
    "# Load sample data from a WAV file\n",
    "sample_rate, data = scipy.io.wavfile.read('ecg.wav')\n",
    "times = np.arange(len(data))/sample_rate\n",
    "\n",
    "# import neurokit2 as nk  # Load the package\n",
    "# ecg = nk.ecg_simulate(duration=10, heart_rate=70)\n",
    "# sample_rate = 70\n",
    "# data = ecg\n",
    "\n",
    "# Apply a 50 Hz low-pass filter to the original data\n",
    "filtered = lowpass(data, 50, sample_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code used to display the result\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3), sharex=True, sharey=True)\n",
    "fig.set_dpi(800)\n",
    "ax1.plot(times, data)\n",
    "ax1.set_title(\"Original Signal\")\n",
    "ax1.margins(0, .1)\n",
    "ax1.grid(alpha=.5, ls='--')\n",
    "ax2.plot(times, filtered)\n",
    "ax2.set_title(\"Low-Pass Filter (50 Hz)\")\n",
    "ax2.grid(alpha=.5, ls='--')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize R-peaks in ECG signal\n",
    "# plot = nk.events_plot(rpeaks['ECG_R_Peaks'], ecg_signal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# High-Pass Filter\n",
    "import numpy as np\n",
    "import scipy.signal\n",
    "import scipy.io.wavfile\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def highpass(data: np.ndarray, cutoff: float, sample_rate: float, poles: int = 5):\n",
    "    sos = scipy.signal.butter(poles, cutoff, 'highpass', fs=sample_rate, output='sos')\n",
    "    filtered_data = scipy.signal.sosfiltfilt(sos, data)\n",
    "    return filtered_data\n",
    "\n",
    "# Load sample data from a WAV file\n",
    "sample_rate, data = scipy.io.wavfile.read('ecg.wav')\n",
    "times = np.arange(len(data))/sample_rate\n",
    "\n",
    "# Apply a 20 Hz high-pass filter to the original data\n",
    "filtered = highpass(data, 20, sample_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code used to display the result\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3), sharex=True, sharey=True)\n",
    "fig.set_dpi(800)\n",
    "ax1.plot(times, data)\n",
    "ax1.set_title(\"Original Signal\")\n",
    "ax1.margins(0, .1)\n",
    "ax1.grid(alpha=.5, ls='--')\n",
    "ax2.plot(times, filtered)\n",
    "ax2.set_title(\"High-Pass Filter (20 Hz)\")\n",
    "ax2.grid(alpha=.5, ls='--')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Band-Pass Filter\n",
    "import numpy as np\n",
    "import scipy.signal\n",
    "import scipy.io.wavfile\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def bandpass(data: np.ndarray, edges: list[float], sample_rate: float, poles: int = 5):\n",
    "    sos = scipy.signal.butter(poles, edges, 'bandpass', fs=sample_rate, output='sos')\n",
    "    filtered_data = scipy.signal.sosfiltfilt(sos, data)\n",
    "    return filtered_data\n",
    "\n",
    "# Load sample data from a WAV file\n",
    "sample_rate, data = scipy.io.wavfile.read('ecg.wav')\n",
    "times = np.arange(len(data))/sample_rate\n",
    "\n",
    "# Apply a 10-50 Hz high-pass filter to the original data\n",
    "filtered = bandpass(data, [10, 50], sample_rate)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code used to display the result\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3), sharex=True, sharey=True)\n",
    "fig.set_dpi(800)\n",
    "ax1.plot(times, data)\n",
    "ax1.set_title(\"Original Signal\")\n",
    "ax1.margins(0, .1)\n",
    "ax1.grid(alpha=.5, ls='--')\n",
    "ax2.plot(times, filtered)\n",
    "ax2.set_title(\"Band-Pass Filter (10-50 Hz)\")\n",
    "ax2.grid(alpha=.5, ls='--')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Low-Pass Cutoff Frequency\n",
    "import numpy as np\n",
    "import scipy.signal\n",
    "import scipy.io.wavfile\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load sample data from a WAV file\n",
    "sample_rate, data = scipy.io.wavfile.read('ecg.wav')\n",
    "times = np.arange(len(data))/sample_rate\n",
    "\n",
    "# Plot the original signal\n",
    "fig, ax = plt.subplots(1,1)\n",
    "fig.set_dpi(800)\n",
    "\n",
    "ax.plot(times, data, '.-', alpha=.5, label=\"original signal\")\n",
    "\n",
    "# Plot the signal low-pass filtered using different cutoffs\n",
    "for cutoff in [10, 20, 30, 50]:\n",
    "    sos = scipy.signal.butter(5, cutoff, 'lowpass', fs=sample_rate, output='sos')\n",
    "    filtered = scipy.signal.sosfiltfilt(sos, data)\n",
    "    ax.plot(times, filtered, label=f\"low-pass {cutoff} Hz\")\n",
    "\n",
    "plt.legend()\n",
    "plt.grid(alpha=.5, ls='--')\n",
    "plt.axis([0.35, 0.5, None, None])\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Use Gustafsson’s Method to Reduce Edge Artifacts](https://swharden.com/blog/2020-09-23-signal-filtering-in-python/)\n",
    "- Artifacts may appear in the smooth signal if the first or last data point differs greatly from their adjacent points. This is because, in an effort to ensure the filtered signal length is the same as the input signal, the input signal is “padded” with data on each side prior to filtering. The default behavior is to pad the data by duplicating the first and last data points, but this causes artifacts in the smoothed signal if the first or last points contain an extreme value. An alternative strategy is Gustafsson’s Method, described in a 1996 paper by Fredrik Gustafsson in which “initial conditions are chosen for the forward and backward passes so that the forward-backward filter gives the same result as the backward-forward filter.” Interestingly, the original publication demonstrates the method by filtering noise out of an ECG recording."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.signal\n",
    "import scipy.io.wavfile\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load sample data from a WAV file\n",
    "sample_rate, data = scipy.io.wavfile.read('ecg.wav')\n",
    "times = np.arange(len(data))/sample_rate\n",
    "\n",
    "# Isolate a small portion of data to inspect\n",
    "segment = data[350:400]\n",
    "\n",
    "# Create a 5-pole low-pass filter with an 80 Hz cutoff\n",
    "b, a = scipy.signal.butter(5, 80, fs=sample_rate)\n",
    "\n",
    "# Apply the filter using the default edge method (padding)\n",
    "filtered_pad = scipy.signal.filtfilt(b, a, segment)\n",
    "\n",
    "# Apply the filter using Gustafsson's method\n",
    "filtered_gust = scipy.signal.filtfilt(b, a, segment, method=\"gust\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the Results\n",
    "fig, ax = plt.subplots(1,1)\n",
    "fig.set_dpi(800)\n",
    "ax.plot(segment, '.-', alpha=.5, label=\"data\")\n",
    "ax.plot(filtered_pad, 'k--', label=\"Default (Padding)\")\n",
    "ax.plot(filtered_gust, 'k', label=\"Gustafsson's Method\")\n",
    "plt.legend()\n",
    "plt.grid(alpha=.5, ls='--')\n",
    "plt.title(\"Padded Data vs. Gustafsson’s Method\")\n",
    "# ax.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Filter Using Convolution](https://swharden.com/blog/2020-09-23-signal-filtering-in-python/#filter-using-convolution)\n",
    "- An alternative strategy to low-pass a signal is to use convolution. In this method you create a kernel (typically a bell-shaped curve) and convolve the kernel with the signal. The wider the window is the smoother the output signal will be. Also, the window must be normalized so its sum is 1 to preserve the amplitude of the input signal. Note that this method exclusively uses NumPy and does not require SciPy.\n",
    "\n",
    "- There are different for handling data at the edges of the signal, but setting mode to valid deletes insufficiently filtered points at the edges to produce an output signal that is fully filtered but slightly shorter than the input signal. See numpy.convolve documentation for additional information.\n",
    "\n",
    "- The kernel shape affects the spectral properties of the filter. Commonly called window functions, these different shapes produce filtered signals with different frequency response characteristics. The Hanning window is preferred for most general purpose signal processing applications. See FftSharp for additional information about the pros and cons of common window functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.io.wavfile\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load sample data from a WAV file\n",
    "sample_rate, data = scipy.io.wavfile.read('ecg.wav')\n",
    "times = np.arange(len(data))/sample_rate\n",
    "\n",
    "# create a Hanning kernel 1/50th of a second wide\n",
    "kernel_width_seconds = 1.0/50\n",
    "kernel_size_points = int(kernel_width_seconds * sample_rate)\n",
    "kernel = np.hanning(kernel_size_points)\n",
    "\n",
    "# normalize the kernel\n",
    "kernel = kernel / kernel.sum()\n",
    "\n",
    "# Create a filtered signal by convolving the kernel with the original data\n",
    "filtered = np.convolve(kernel, data, mode='valid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the result\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 3))\n",
    "fig.set_dpi(800)\n",
    "ax1.plot(np.arange(len(kernel))/sample_rate, kernel, '.-')\n",
    "ax1.set_title(\"Kernel (1/50 sec wide)\")\n",
    "ax1.grid(alpha=.5, ls='--')\n",
    "\n",
    "ax2.plot(np.arange(len(data))/sample_rate, data)\n",
    "ax2.set_title(\"Original Signal\")\n",
    "ax2.margins(0, .1)\n",
    "ax2.grid(alpha=.5, ls='--')\n",
    "\n",
    "ax3.plot(np.arange(len(filtered))/sample_rate, filtered)\n",
    "ax3.set_title(\"Convolved Signal\")\n",
    "ax3.margins(0, .1)\n",
    "ax3.grid(alpha=.5, ls='--')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
